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We show how to employ thermal lattice gas models to describe non-equilibrium phenomena. This is 
achieved by relaxing the restrictions of the usual micro-canonical ensemble for these models via the 
introduction of thermal "demons" in the style of Creutz. Within the Lattice Boltzmann approxima- 
tion, we then derive general expressions for the usual transport coefficients of such models, in terms 
of the derivatives of their equilibrium distribution functions. To illustrate potential applications, 
we choose a model obeying Maxwell-Boltzmann statistics, and simulate Rayleigh-Benard convection 
with a forcing term and a temperature gradient, both of which are continuously variable. 
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I. INTRODUCTION 

Frisch, Hasslacher and Pomeau (FHP) pioneered 
the use of Lattice Gas Automata (LGA) to simulate the 
Navier-Stokes (NS) fluid. The motion of fictitious parti- 
cles on an underlying hexagonal lattice, subject to care- 
fully chosen rules for collisions and propagation, gives 
rise to the NS equations in the continuum limit. Since 
that time, the LGA model and its derivative, the Lattice 
Boltzmann (LB) model have attracted considerable 
attention because of their potential application to the 
simulation of complex fluid systems, in particular, sys- 
tems with inter-particle interactions which model phase 
transitions and the dynamics of interfaces. Although the 
earliest such models achieved spatial variation of the 
order parameter only by ignoring even semi-detailed bal- 
ance in the description of possible collisions, more recent 
models are free from such inconsistencies 

However, like the original FHP model, all these models 
are intrinsically "athermal" Q , and so are unable to sim- 
ulate phenomena where the temperature is an important 
variable. Only recently have thermal LGA models been 
constructed, such that the thermodynamics of fluids can 
be studied ||7[-pT|, and, typically, they are defined within a 
micro-canonical ensemble. To satisfy strict conservation 
rules for particle number, momentum and energy, they 
must contain several species of particles with different 
energies, interacting via rules which are often very com- 
plex. Reference Q, for example, is an extension of the 
FHP model to 4 species of particles with carefully cho- 
sen momenta and energies, and in references the 
authors have to introduce unequal masses for particles 
moving in different channels, in order to ensure a suffi- 
cient number of allowed collisions. Thermal LB models 
are in principle less complex because they are expressed 
in terms of the distribution functions for the fictitious 
particles, but have not been conspicuously successful 

A more serious problem for the study of thermal 
phenomena, particularly those of a non-equilibrium na- 
ture (heat conduction) or involving instabilities (convec- 



tion), arises because both LGA and LB thermal models 
treat the temperature as an externally defined param- 
eter. However, when the temperature is itself a spa- 
tially varying parameter, it is necessary to have infor- 
mation about its variation in order to implement the col- 
lision rules (LGA) or the relaxation process (LB). Con- 
sequently, when this spatial variation is itself the object 
of study, it is not possible without modification p^ ] to 
apply existing models. 

Our purpose in this paper is to introduce a class 
of thermal models which permit the study of non- 
equilibrium phenomena. As a by-product, these models 
permit a much wider variety of collision rules for LGA 
models, and are also readily adapted to the LB approach. 
The key feature is a novel idea drawn from the Monte- 
Carlo literature |p^|l6|. In the language of the Ising 
model, each site is associated with a local thermal reser- 
voir or "demon" , which interchanges energy with the par- 
ticles on the site. The local thermodynamic temperature, 
in the low velocity limit, is then proportional to the lo- 
cal demon energy. In the LGA or LB models, each de- 
mon is associated with one node of the hexagonal lattice, 
and constitutes a mechanism for monitoring or control- 
ling the local temperature, so that the modeling of non- 
equilibrium phenomena becomes possible. 

In the following section, we describe the equilibrium 
properties of such thermal lattice gas models. Then, 
in Section 
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we derive expressions for the transport 
coefficients of Lattice Boltzmann models by means of a 
Chapman-Enskog expansion. In Section |rv| we illustrate 
ideas using a particular model, and in Section^ as a spe- 
cific example, use this model to simulate two-dimensional 
Rayleigh-Benard convection and to draw some prelim- 
inary conclusions concerning the feasibility of our ap- 
proach. 



II. THERMAL LATTICE GAS MODELS 
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A. Generalities 

In analogy with the basic LGA model Q| , we define a 
class of thermal models in which several species of fic- 
titious particles move on a two-dimensional hexagonal 
lattice. Some of the particles (rest particles) remain at 
rest. During each time step, the particles interact ("col- 
lide" like billiard balls) and then propagate ballistically 
at constant speed from one site of the lattice to a neigh- 
bouring site. The collision step rigorously conserves the 
number of particles and their total momentum. Energy 
is also conserved, but in a particular manner, described 
below. 

Energy is defined as the sum of the kinetic and poten- 
tial energies of the individual fictitious particles, and it 
is essential that there be at least two species of particles 
with distinct energies for there to be well-defined thermal 
properties. In the literature, such thermal models often 
employ only kinetic energies, so that within the micro- 
canonical ensemble, particles with different speeds and 
or masses Jsj-p^ are required to satisfy the conservation 
laws. (This restriction to kinetic energies and to colli- 
sion processes which strictly conserve energy also has the 
consequence jl^ that such models posses zero bulk vis- 
cosity.) However, the inclusion of potential energy terms 
is straightforward, |l^J^, although seldom employed for 
other than pedagogical purposes. 

Conservation of energy is enforced with the aid of the 
so-called demons, one at each site. The advantage of 
this procedure is that it is not necessary to choose values 
of masses, speeds and/or energies to permit a sufficient 
number of non-trivial collisions. Instead, a demon acts 
as a kind of energy reservoir, permitting collisions in a 
manner which satisfies detailed balance. Over the course 
of time, there will be a distribution function for the local 
demon energy, which has the form P{Ed) ^ exp —PEjj, 
where /3 = 1/T is the inverse of the local temperature. 
Since Ejj is a continuously variable quantity, its average 
value < Ed > is equal to the local temperature. The 
idea is borrowed from Creutz p^ , and was exploited in 
a series of papers applying Ising-style lattice gas models 
to non-equilibrium interface problems [l^ . 

At each site, a collision process may proceed either 
if it produces surplus (or zero) energy or if the local 
demon can provide the requisite energy deficit. In the 
former case, the demon absorbs the surplus. Demons 
are required always to have positive energy, and as a re- 
sult serve also to regulate the occurrence of collision pro- 
cesses. As demonstrated by Creutz and by Jorgenson et 
al 1 18 1, the demons thus act as thermometers, measuring 
the local temperature by virtue of their own statistically 
averaged energy. It is also possible to use the demons 
to control the local temperature: if this is done at the 
boundaries of a sample, then, for example, a tempera- 
ture gradient can be set up across the sample, permit- 
ting a measurement of thermal conductivity pof or the 
establishing of Rayleigh-Benard convection. 



A useful way to understand the role of the demons 
is to consider them as a species of rest particle. Rest 
particles of the conventional kind can be created or anni- 
hilated in collision processes so as to satisfy the conser- 
vation laws. Demons are neither created nor annihilated, 
but act to satisfy the conservation of energy. Conven- 
tional rest-particles (usually) have only one energy level, 
but demons have many such levels - although in simple 
situations one might ima gin e demons with only a small 
number of distinct levels 



B. Statistical Equilibrium 

As a first step in deriving expressions for the transport 
coefficients we define some generic notation Sites in 
the hexagonal lattice are labeled by an index i and by 
a site vector R^. The vectors radiating from a site i to 
its (not-necessarily) nearest neighbours define different 
possible directions for the motion of fictitious particles 
at that site. For particles of species /, we have c^,/, a — 
1,2,- • • , 67, where hi is the number of distinct neighbours 
at distance c/. This distance also defines the "speed" of 
the species /, so that the kinetic energy of each fictitious 
moving particle is \(?i- In similar fashion, each particle of 
species / has potential energy ej, so that its total energy 
is El = ^cj + €i. Any rest particle has total energy zero. 

The fictitious particles occupy the available states 
{i, a, 1} according to particular statistics, and the de- 
tailed balance property of the collision rules guarantees 
that the system possesses a state of statistical equilib- 
rium. Furthermore, the existence of an equilibrium dis- 
tribution function is an essential condition for the regain- 
ing of the Navier-Stokes equations in the continuum limit 
The usual choice of statistics is Fermi-Dirac, since, 
historically, it was convenient for coding purposes to rep- 
resent occupation numbers as binary variables. However, 
this choice is not essential. In general we write the oc- 
cupation number of the state {i,a, /} as rii^aj, and the 
ensemble-averaged distribution function for this state as 
fi,a,i =<ni^a,i >■ It is convenient to represent rest par- 
ticles in the same way by rii^ and fifl. 

In an equilibrium state for which there is no net mo- 
tion of the fiuid, the functions / have the values /j', so 
that the probability of finding a certain configuration can 
then be written as a product of the continuous variables 
/^^ and P{Eo): the latter variable being the probability 
of finding a demon with energy Eo. It is now possible to 
define thermodynamic variables in terms of the equilib- 
rium distributions. 

The most important of these are the density p and the 
pressure, P. The density is just 



PI 
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where M is the number of rest-particle states per site. 
In principle, the pressure should be derived from the free 
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energy. However, as is well known in the lattice gas liter- 
ature , the quantity which plays the role of pressure 
in the hydrodynamic equations is 



P 



2 ^ 



(2) 



sometimes known as the "kinetic pressure" . It is useful 



to define the analogous partial pressures Pj 



1 ^2 



I Pi, 



that P = We will also require the energy den- 

sity U — J2i Eipi and the isothermal and adiabatic 
speeds of sound, which are respectively 



4 = \ ^Pic]Ip 



(3) 



and 



interact with each other ("collide") following predeter- 
mined rules 1^. Energy conservation is ensured by the 
local demon, as described previously, and detailed bal- 
ance is rigorously observed. Conservation of particles, 
momentum and energy on each site can be written ex- 
plicitly as 



fi,a,i{t) + fifi{t) = constant 



a, I 



Cafi,a.i{t) = constant 



fi,a,i{t)Ei + Eo^it) = constant 



(10) 
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(12) 
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The derivatives of the equilibrium distribution func- 
tions are defined as ^ T (^^^ where T is the 

temperature and p, is the chemical potential | |2l[ |. 

When there exists a net local flow velocity u it is still 
possible to define equilibrium distribution functions. In 
the low velocity limit u ^ c, they can be expanded in 
powers of u. Following reference , but with a slightly 
generalised notation, we obtain to first order p2| 



fZ{n) = fr + qpf^c,J-n 



(5) 



The constant q is determined by requiring that the ex- 
pansions be consistent with the total momentum defined 
as 



■J 



We obtain 



g = 2/V6,/^? 



and it is convenient to define qi = qfj so that 
/aM") = fr + P<liCaJ ■ u 



(6) 

(7) 

(8) 

(9) 



and i ^bjqjcj = 1. Note that to first order in u, there 
are no corrections to /q' or to P{Eo)- 



where t is the (discrete) time. E^^i is the demon energy. 

After each collision step, there is a propagation step, in 
which each moving particle moves one lattice constant in 
the direction of its velocity, while the rest particles and 
the demon remain unchanged. Thus the kinetic equa- 
tions for the particle and demon variables, including both 
collisions and propagation, are: 



/j,M(*+l)-Aa,/(i)=f^/, 

hfl{t + 1) - hflit) = 



1, 



(13) 



(14) 



where the index j is defined by Rj — + c^. ^ is the 
ensemble averaged collision operator. Strictly speaking, 
ED^i should also be replaced by its ensemble average, 

<ED,i>. 

The lattice Boltzmann approximation to the collision 
matrix directly employs the existence of equilibrium dis- 
tribution functions. Any distribution function which 
deviates from its equilibrium value can be written as 
/ /'^'^ + 5j where g is hopefully small. If we assume 
that there exists a single characteristic relaxation time 
T [p3| , then in terms of g, we may approximate the rate 
equations as 

f],a,lit+l) - fi.ajit) ^ ~gi,a.l/T, 0=1,- ••,5/ (15) 
f^At+^)-koit)^-9^,o/T (16) 

These expressions will be employed to simplify the anal- 
ysis in the subsequent parts of the paper. 



C. Time Evolution 

As in a traditional LGA , time evolution proceeds in 
two distinct steps. First, particles at a given site, (both 
rest particles and moving particles in any energy level). 



III. TRANSPORT COEFFICIENTS 

The basic continuum equations for a thermal fluid 
could now be obtained via a Chapman-Enskog expansion. 
However, since such treatments are available elsewhere 
|lOt , we will focus only on the derivation of expressions 
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for the transport coefficients. We will follow the pro- 
cedure of reference |2|1. The first step is to replace the 
rate equations (|l^), ( |16D by their continuum versions. In 
the limit where spatial and temporal changes are small 
and/or slow, the left-hand sides of these equations yield 
derivatives which can be replaced by dt — edt^ +e^dt2, 
Vr = eVri, where e is a small parameter. The distri- 
bution functions can also be expanded in terms of e. To 
order e^, this gives: 



and 



Ja,I - Ja,I + <^Ja,I + /a,/' 



f - f(o) 4_<:y-(i) _U^2 .(2) 
/o,7 - Jo,/ + e/oj + e hj 



where the zeroth order terms in the expansion are just 
the equilibrium distributions, and the higher order terms 
are just g: /^"^ = /(°) and g = e/^^' +t^f^'^\ Necessarily, 
the conservation of mass, momentum and energy require 
that 



and 



E 



{a) 



(17) 
(18) 

(19) 



and 



(26) 



The expressions for the shear viscosity ly, the bulk vis- 
cosity C and the thermal conduction coefhcient A arise 
from the second order terms in e. After some reduction, 
the corresponding equations become 



^*2(p) + iE<Ki/l°i = o 



(27) 
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and 



a, I 



a, I 



(28) 



where the terms preceded by the factor i are second order 
in the Taylor expansion of the derivatives. These terms 
can be reduced by makinguse of the equations ( |l7|) and 
( pO| ) . Although equation ( p7| ) reduces to the simple result 
9(2 (p) = 0, equation ( ^ is more interesting. We simphfy 
its second term to read X]a/^a,i ■ ViI?^^)caj/^°], and 
then substitute for f^^j from ( pO| ) to obtain 



9,2P(U) = (r - i) E Ca./ • Vii?i>a,l/S 



(29) 



where a = 1,2. 

Substituting into equations (|l^)-(16), and equating 
terms for each order of e, we obtain a hierarchy of Boltz- 
mann equations. The first order equations are: 
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(20) 






(21) 






(22) 



Using ( |l7|) and (|l^) we then obtain 
9ti(p)+Vi-(pu) =0 



and 



dti{pu) + \/iP = 



(23) 



(24) 



where P is the kinetic pressure as defined earlier. Simi- 
larly, since these relations are true for any values of the 
distribution functions /^''(u = 0), there must also be 
analogous relations for the partial densities pi and par- 
tial pressures P/, namely 



aa(p/) + — Vi-(pu) = 
P 



(25) 



The final step is to substitute (^, and use ( psj ) and ( |2q ) 
to write 

dt2{pu) = (r - i) ^(X/V^u + 17 VV • pu) (30) 
where 



Xr 



rbiqicj 



and 



Yi = -^{biqicj - pic]/p) 

Combining with (|2^), this becomes 

dt{pn) = - VP + (r - i)(z^V2pu + (VV • pu) (31) 

where the shear viscosity u is 

1 



and the bulk viscosity Q is 



C=(c^-4)(r-i) 
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using the expressions for the sound speeds Ct and Cg 
given earUer. Note that the bulk viscosity vanishes for 
any model having only one value of the speeds c/: this 
includes any FHPl model, ||], as a special case. 

The procedure to obtain the thermal conduction coef- 
ficient, A, is very similar to that described, for example, 
in Huang |^^. A key feature is the necessity to subtract 
out of the energy current that part which depends on 
the net flow of the fluid. This is the origin of the "sub- 
tracted current" described by Ernst QJl^ , but neglected 
by Chen et al ||T^ . We account for this effect by replacing 
the energies Ei by Ej such that 



(32) 



The 



where /3 = i and fj is given hy fj — — 

high temperature behaviour of A is dominated by the yj- 
prefactor. 

For the particular case of the Fermi-Dirac distribution, 
it is easy to show that this is precisely the expression 
given by Ernst Q- We write /f = ~{fi-Ei)f'^ and 
therefore obtain 

^ = - J)^ E bif^cjEjiEi - n)/2 



2'T' 
1)1. 



J2bif^cjEj/2 



(40) 



where the last step comes from the identity (k3 



or 



I 



bifj^c^jEi = 



(33) 



Writing Ei = Ej — Es, we obtain 

Es = ^ bjqiEicj/ ^ biqjcj 
I I 

= J2bjf^Ejcj/J2bjf^cj 



(34) 



The first order equation describing the time evolution 
of the local energy density U ~ J^a i fa'^ii^)^! then 



dti[U + ED] = 



(35) 



which gives no new information. However, after using 
equations (19)-(^0|), the second order equation becomes 



1, 



dt2 [U + En] = {r - EiCaj ■ Vi^« /;«„(u) (36) 

a, I 



Using the explicit form of /^^^(u), equation (^, and the 
definition of E, equation (|33|), we therefore obtain 

dt[U + ED] = {T-l)Y,EicjW^Pi (37) 



It remains to express the gradient explicitly in terms of 
the thermal gradient. In the absence of a net flow of 
particles, the chemical potential ^ is not a function of 

position, and so V^/j = ff\7^T, where /f = , 
and the coefficient of thermal conduction is 



r-hY.bifIc'jEj/2 



It is convenient to write this also in the form 
1, 1 



A=(t--) 



E bif?cjEj/2 



(38) 



(39) 



IV. A MODEL 

To illustrate the ideas of the previous sections, and, 
in particular, to illustrate the use of statistics other than 
Fermi-Dirac, we have previously described a model where 
the statistics are quite unconventional [ p5[ . In the present 
paper, we choose to employ Maxwell-Boltzmann statis- 
tics for a simple model with 3 energy levels. The lowest 
level, with energy Eq = 0, is M-fold degenerate, so that 
there are at most M fictitious particles with zero velocity 
("rest particles"). The other two levels both correspond 
to the same speed c, with the same 6-fold symmetry as 
in the FHP models, but their energies are respectively 
Ea and Eb = Ea + A. In our numerical calculations, we 
take Ea — 0.62 and Eb = 1-80 in units of c^. 

Assuming Maxwell-Boltzmann statistics, we obtain 



/o(u) 



//,a(u) 



(l + qpu-Ca), I^A,B 



where q = pe~^''/3c2(e"^^-^ + e-/^^^), and pe^'^'" = 
[M+6{e^'^^^+e^'^^^)]. This simple form for the velocity 
expansion results because fj = fj"^. 

Because there is only one speed in the model, the ex- 
pressions for the sound velocities and for the shear and 
bulk viscosities become particularly simple. We obtain 



C'Y' 



and 



1 2^a,bPi 



—c 
2 



3c' 



^-I3Ea _|_ (,-I3Eb 



M -f 6(e-''^'^ +e-P^'^) 



(41) 



V = [t 
= {r 



>V4; 

2' 2 p 
l.c2 



M 



2' 2 M + Gie-f^^^ 



(42) 



The temperature-independence of v is an artifact of the 
model, but the temperature dependence of Q is typical of 
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the viscosity of a liquid. For a typical value of the density, 
C is displayed in Figure 1 as a function of temperature. 

The model also leads to a simple expression for the 
thermal conduction coefficient A: 



A = (r- 



1 ScW 

r T2 



neq req 
J A J B 



req . req 
J A + Jb 



For the purpose of illustration, it is convenient to relate 
this to the thermal diffusion coefficient, defined as Dt = 
X/ pCp, which is also plotted in Figure 1. Although the 
temperature dependence of A is dominated by the ^ 
prefactor at high temperatures, T > A, this behaviour is 
exactly compensated by the temperature dependence of 
the specific heat. 
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FIG. 1. Temperature variation of the transport coeffi- 
cients for the model described in the text. Solid line: bulk 
viscosity. Dashed line: shear viscosity. Dotted line: thermal 
diffusion. Dash-dotted line: sound attenuation. All coeffi- 
cients are in units of r — i. 

For completeness, Figure |l| also shows the temperature 
variation of the sound attenuation coefficient F, defined 
in the standard way as F = + C + (7 ~ 1)^t)- 



V. RAYLEIGH BENARD CONVECTION 

To illustrate the feasibility of our approach, we car- 
ried out simulations of Rayleigh-Benard convection. We 
considered a two dimensional cell with horizontal length 
L and vertical height H . Periodic boundary conditions 
were imposed in the horizontal direction, with two rigid 
walls at the top and the bottom of the cell. The de- 
mon energies on the upper and lower walls were fixed 
at values Ti and T2, respectively. A uniform force was 
implemented by changing the vertical momentum of the 



particles at a constant rate, while keeping horizontal mo- 
mentum unchanged. Thus, both the temperature differ- 
ence 5T and the force could be tuned continuously. 

Our initial study was of a system with L — 400 lattice 
units and H — 100 x VS units, so that the aspect ratio 
L : H was near 2 : 1, with a particle density of 3.6 per 
site. The energy levels were Ea = 0.62, Eb — 1-8 in 
units of c^, and the temperatures at the lower and upper 
boundaries were Ti = 4.8, T2 = 0.3 in the same units. 
We chose the relaxation time r to be 1 unit, and averaged 
velocity fields over 20 x 20 regions in order to evaluate 
the local distribution functions. Initial conditions were a 
uniform density and a uniform temperature gradient, and 
with the local velocity everywhere zero, so that we were 
able to observe the onset and subsequent evolution of the 
convective instability. With these parameters, the evolu- 
tion of the system was sufficiently slow that "snapshots" 
of the velocity field could be obtained by time averages 
over only 50 time-steps. 

Figure ^ shows the evolution of the temperature distri- 
bution, and of the corresponding velocity fields for these 
parameters. At first four convection rolls are clearly seen, 
but subsequently these collapse into two. Evidently, the 
model has captured the essential features of the physical 
phenomenon. Systematic analysis of the data as a func- 
tion of the system parameters will be the subject of a 
subsequent publication, but certain preliminary remarks 
are appropriate here. 

According to the classic linear stability analysis, jp6[, 
convection occurs when the Rayleigh number TZ ||27| ex- 
ceeds a critical value of order 10^. Even for situations far 
from the linear regime, the value of 7?. is a good indica- 
tor of the stability of the convective phenomenon. Thus, 
since TZ is of order 10^ for the simulation shown in the fig- 
ure, the convection rolls should be extremely stable, as is 
indeed the case. Indeed, our results are strikingly similar 
to those of a previous detailed study |^ , obtained with 
a modified LB model which represented the temperature 
as a "convected passive scalar field" . 

An interesting feature of our results is the observa- 
tion of four rolls before the final stable state with two 
rolls is reached. Linear analysis predicts two rolls (since 
the wavenumber of the instability should be around ir/H 
where H is the height of the cell) . Our tentative explana- 
tion is that the velocity field is first established near the 
lower (hot) boundary, so that the effective height of the 
cell is considerably smaller than H. We therefore expect 
the wavenumber of the instability to be larger than vr/iJ, 
and the periodic boundary conditions select ^ 2tt/ H . (In 
principle, with different geometry, we might expect to see 
yet higher initial wavenumbers.) The study of this tran- 
sient regime, and the subsequent stabilisation of the two 
rolls, will form part of our ongoing investigations. 
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FIG. 2. Contour plots showing the time evolution of the 
temperature distribution (left panels) and of the velocity field 
(right panels). From top to bottom, data is shown at times 
20,000, 40,000, 60,000, 80,000 and 100,000 units respectively. 
One time unit is defined as one update of every site in the 
lattice. 

However, in general terms, this preliminary illustration 
of Rayleigh-Benard convection demonstrates the feasibil- 
ity of our Lattice Boltzmann method for non-equilibrium 
thermal lattice-gases. We intend to apply our approach 
to a variety of phenomena for which the statistical me- 
chanics of the gas are critically important. In particular, 
we plan to include interactions between fictitious par- 
ticles so as to simulate systems with first order phase 
transitions and the dynamics of interfaces between their 
associated phases. 

We thank Hong Guo and Martin Grant for many useful 
discussions. 
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